####################################################
# Generalized DiD: Crime Concerns Other Countries
####################################################

rm(list=ls())

library(Hmisc)
library(readstata13)
library(foreign)
library(lmtest)
library(sandwich)
library(stargazer)
library(msm)
library(ggplot2)
set.seed(1111)

#sink("14_generalized_did_concerns_other_countries.txt")

# functions
vcovCluster <- function(
  model,
  cluster
)
{
  require(sandwich)
  require(lmtest)
  if(nrow(model.matrix(model))!=length(cluster)){
    stop("check your data: cluster variable has different N than model")
  }
  M <- length(unique(cluster))
  N <- length(cluster)           
  K <- model$rank   
  if(M<50){
    warning("Fewer than 50 clusters, variances may be unreliable (could try block bootstrap instead).")
  }
  dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
  uj  <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum));
  rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
  return(rcse.cov)
}

round_df <- function(df, digits) {
  nums <- vapply(df, is.numeric, FUN.VALUE = logical(1))
  
  df[,nums] <- round(df[,nums], digits = digits)
  
  (df)
}

############################
# Read and prepare data 
############################

# load individual data
load("final_individual_2023april6.RData")
names(d)

# binary treatment
table(d$t_stagg_venezuela, exclude = NULL)
table(d$t_stagg_colombia, exclude = NULL)
table(d$t_stagg_bolivia, exclude = NULL)
table(d$t_stagg_argentina, exclude = NULL)
table(d$t_stagg_ecuador, exclude = NULL)

############################################
# Concerns: binary treatment covariates
############################################

# Venezuela
reg1 = lm(d$crime_concerns_s ~ d$t_stagg_venezuela + d$gender_respondent + d$age_bins_respondent + d$education_respondent + as.factor(d$county) + as.factor(d$survey))
r1 = round(coeftest(reg1, vcov = vcovCluster(reg1, cluster = d$county)),3)

# Colombia
reg2 = lm(d$crime_concerns_s ~ d$t_stagg_colombia + d$gender_respondent + d$age_bins_respondent + d$education_respondent + as.factor(d$county) + as.factor(d$survey))
r2 = round(coeftest(reg2, vcov = vcovCluster(reg2, cluster = d$county)),3)

# Bolivia
reg3 = lm(d$crime_concerns_s ~ d$t_stagg_bolivia + d$gender_respondent + d$age_bins_respondent + d$education_respondent + as.factor(d$county) + as.factor(d$survey))
r3 = round(coeftest(reg3, vcov = vcovCluster(reg3, cluster = d$county)),3)

# Argentina
reg4 = lm(d$crime_concerns_s ~ d$t_stagg_argentina + d$gender_respondent + d$age_bins_respondent + d$education_respondent + as.factor(d$county) + as.factor(d$survey))
r4 = round(coeftest(reg4, vcov = vcovCluster(reg4, cluster = d$county)),3)

# Ecuador
reg5 = lm(d$crime_concerns_s ~ d$t_stagg_ecuador + d$gender_respondent + d$age_bins_respondent + d$education_respondent + as.factor(d$county) + as.factor(d$survey))
r5 = round(coeftest(reg5, vcov = vcovCluster(reg5, cluster = d$county)),3)

############################################
# Concerns: binary treatment no covariates
############################################

# Venezuela
reg6 = lm(d$crime_concerns_s ~ d$t_stagg_venezuela + as.factor(d$county) + as.factor(d$survey))
r6 = round(coeftest(reg6, vcov = vcovCluster(reg6, cluster = d$county)),3)

# Colombia
reg7 = lm(d$crime_concerns_s ~ d$t_stagg_colombia + d$gender_respondent + as.factor(d$county) + as.factor(d$survey))
r7 = round(coeftest(reg7, vcov = vcovCluster(reg7, cluster = d$county)),3)

# Bolivia
reg8 = lm(d$crime_concerns_s ~ d$t_stagg_bolivia + as.factor(d$county) + as.factor(d$survey))
r8 = round(coeftest(reg8, vcov = vcovCluster(reg8, cluster = d$county)),3)

# Argentina
reg9 = lm(d$crime_concerns_s ~ d$t_stagg_argentina + as.factor(d$county) + as.factor(d$survey))
r9 = round(coeftest(reg9, vcov = vcovCluster(reg9, cluster = d$county)),3)

# Ecuador
reg10 = lm(d$crime_concerns_s ~ d$t_stagg_ecuador + as.factor(d$county) + as.factor(d$survey))
r10 = round(coeftest(reg10, vcov = vcovCluster(reg10, cluster = d$county)),3)

###############################
# Tables
###############################

# Table a12
stargazer(r6,r1,r7,r2,r8,r3,r9,r4,r10,r5,
          digits = 3,
          title="",
          align=TRUE, 
          omit.stat=c("LL","ser","f","rsq","adj.rsq"), 
          no.space=TRUE,  
          multicolumn = TRUE,
          table.placement = "H",
          column.separate = c(2,2),
          covariate.labels=c("Venezuela","Colombia","Bolivia","Argentina","Ecuador"),
          dep.var.caption  = "Crime Concerns" ,
          dep.var.labels.include = FALSE,
          star.cutoffs = c(0.05,0.01,0.001),
          omit = c("Constant","gender_respondent","age_bins_respondent","education_respondent","county","survey"),
          add.lines = list(c("Covariates","No","Yes","No","Yes","No","Yes","No","Yes","No","Yes"),c("Observations","13,551","13,551","13,551","13,551","13,551","13,551","13,551","13,551","13,551","13,551")),
          out="tableA12.txt")

sink()
